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I. 


INTRODUCTION 


This  work  represents  a continuation  of  an  effort  begun 
considerably  earlier  by  Mullinix  and  Smith  [1].  That  report  dealt 
with  surface  flaws  in  hollow  cylinders  which  were  subjected  to  beam 
bending  action.  The  present  work  also  deals  with  surface  flaws  in 
cylinders,  but  here  the  cylinders  are  subjected  to  two  additional 
loading  cases:  uniaxial  extension  and  internal  pressurization. 
Additionally  for  the  latter  loading,  two  surface  flaw  orientations  are 
considered:  the  flaw  oriented  transverse  to  the  cylinder  axis  as  in 

[1]  and  the  flaw  oriented  in  the  longitudinal  direction. 

As  noted  in  [1],  the  surface  flaw  is  of  considerable  technological 
interest.  It  lias  caused  a number  of  well-documented  service  failures. 
These  failures  can  be  catastrophic  and,  depending  on  the  severity  of 
the  flaw,  can  occur  at  a fairly  low  stress  level.  This  situation  is 
aggravated  by  the  ever-present  demand  for  lighter  weight  structures, 
particularly  in  missile,  rocket,  and  aircraft  applications,  Because 
of  its  importance,  the  surface  flaw  was  the  subject  of  a special 
symposium  at  the  1972  annual  meeting  of  the  American  Society  of 
Mechanical  Engineers  (ASME)  [2].  A collection  of  surface  flaw  solutions 
was  included  in  the  fracture  handbook  by  Smith  and  Mullinix  [3]  as  well 
as  in  the  report  by  Paramerter  [4].  Still,  insufficient  surface  flaw 
solutions  exist,  and  experimental  efforts  to  obtain  solutions  for 
specific  cases  continue;  as  in  the  recent  article  by  Smith,  et  al.  [5], 
for  example. 

Surface  flaws  normally  occur  in  a semi-elliptical  shape.  In  order 
to  cut  flaws  in  the  cylinders  with  circular  saws,  the  flaws  in  this 
report  are  represented  as  part  of  a circle  penetrating  various  depths 
into  the  walls  of  the  cylinders.  There  appear  to  be  no  significant 
errors  involved  in  using  these  particular  flaw  geometries. 


II.  ANALYTICAL  CONSIDERATIONS 


The  equations  of  the  stresses  near  the  edge  of  a surface 
flaw  are  now  well  known  [6,7,8].  In  fact,  in  a plane  perpendicular 
to  the  crackfront,  the  stress  distribution  near  the  part-circular 
crack,  corresponding  to  the  opening  mode  of  deformation,  can  be 
taken  as 
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where  is  the  stress  intensity  factor  and  the  coordinates  r and  iJj 

are  as  shown  in  Figure  1.  Thus,  the  effect  on  the  stress  field  of  the 
crack  border  curvature  as  well  as  the  location  and  shape  of  other 
boundaries  is  reflected  in  the  magnitude  of  K , the  stress  intensity 

factor.  The  stresses  given  by  Equation  (1)  are  valid  only  very  near 
the  crack  tip  where  it  is  assumed  that  the  stresses  singular  in  r are 
much  larger  than  other  terms  which  are  regular  in  r.  Experimental 
measurements,  however,  must  be  made  at  a finite  distance  from  the 
crack  tip  where  other  terms  may  make  a significant  contribution  to 
the  total  stress  field.  These  other  terms  must  somehow  be  accounted 
for  if  Kj  is  to  be  determined.  A simple  way  of  approximately  accounting 

for  the  regular  terms  was  provided  by  Irwin  [9].  Irwin  assumed  a 

uniform  stress  field,  a , at  the  crack  tip  parallel  to  the  crack 

on  r 

plane.  With  this  superposition,  the  form  of  the  local  stress  field 
becomes 


Ki  ± 

O = TJo  COS  9 

" (2irr) 1/2  2 

K 


1 sin  ^ 

2 2 


y <2«) 


lTn  cos  ! (2  + sln  f sin  ¥ 


)- 

) 


on 


ny 


(2nr) 


. ♦ 4)  34> 

jj2  sln  2 cos  2 cos  2 


(2) 


While  o has  no  influence  upon  the  singular  stress  field  itself, 

it  does  alter  the  isochromatic  fringe  pattern,  which  is  proportional 
to  the  maximum  in-plane  stearing  stress.  In  the  use  of  photoelasticity, 
the  most  readily  obtained  quantity  is  the  maximum  shearing  stress, 

Tmax’  From  this  stress  field,  the  maximum  shearing  stress  in  the  plane 

perpendicular  to  the  crack  front,  y-n,  can  be  determined  from 
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Photoelast ically,  the  maximum 
stress  optic  law 
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shearing  stress  is  determined  from  the 


(5) 


where  t is  the  thickness  of  the  specimen  measured  parallel  to  the  direc- 
tion of  the  light,  N is  the  isochromatic  fringe  order,  and  f is  the 
photoelastic  fringe  constant  for  the  material. 


It  is  convenient  to  measure  t along  the  line  = it / 2 where  the 

max 

maximum  shearing  stress  is  known  to  be  large.  Simplifying  Equation  (4) 
at  ip  = n/2  results  in 
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By  measuring  t at  two  different  distances,  r,  from  the  crack  tip,  one 

can  obtain  two  equations  which  can  be  solved  for  the  unknown  parameters 

Kt  and  a . More  than  two  measurement*,  result  in  a statistically  improved 
I on 

result  for  KT  and  a . Only  the  stress  intensity  factor  KT  is  used  in 
I on  I 

fracture  criteria. 


III.  EXPERIMENTAL  CONSIDERATIONS 

The  surface  flaw  problem  is  inherently  three-dimensional  in 
nature  and  requires  a three-dimensional  photoelastic  capability.  Appli- 
cation of  three-dimensional  photoelasticity  to  fracture  mechanics  for 
estimating  is  now  widely  accepted  [5]. 

In  the  present  test  cases,  the  material  used  was  Hysol  CP5-4290 
cast  by  the  Hysol  Corporation,  Olean,  New  York.  The  cylinders  are  nomi- 
nally 5.875  in.  in  outside  diameter  with  a 0.75-in.  w^li  thickness.  Two 
specimen  configurations  were  used:  one  with  the  surface  I law  transverse 
to  the  cylinder  axis  and  one  with  the  flaw  oriented  ii  the  longitudinal 
direction  as  shown  in  Figure  2.  The  flaws  were  machined  with  a circular 
saw  blade  0.006-in.  thick.  Blade  radii  of  0.875  and  1.500  in.  were 
used  to  produce  flaws  of  two  different  sizes. 

As  indicated  in  Figure  2,  a total  of  three  series  of  tests  were 
conducted:  a series  for  the  transversely-flawed  cylinder  loaded  in 

uniform  uniaxial  tension,  a series  for  the  transversely-flawed  cylinder 
loaded  with  internal  pressure,  and  a series  for  the  longitudinally- 
flawed  cylinder  loaded  with  internal  pressure.  The  uniform  tension  load 
was  applied  simply  by  hanging  dead  weights  on  the  cylinder.  The  inter- 
nal pressure  load  was  supplied  by  compressed  air  passing  through  a 
regulator.  The  pressure  was  measured  by  a Mercury  manometer.  Figure  3 
is  a schematic  of  the  setup  for  the  internal  pressure  loading. 

After  the  surface  flaw  was  machined,  the  cylinder  was  annealed  by 
thermal  soaking  at  280°F  for  6 hours  followed  by  cooling  at  the  rate  of 
l°F/hr.  The  stress  freezing  of  the  loaded  models  was  accomplished 
using  the  same  thermal  cycle  as  for  annealing  except  under  loading  con- 
ditions. After  the  stress  freezing  cycle,  slices  perpendicular  to  the 
crack  border  were  removed  from  the  model  by  means  of  a band  saw.  Figure 
4 shows  a few  of  the  various  angles  at  which  slices  were  taken.  The 
number  of  slices  varied  from  test  to  test  depending  upon  the  flaw  size. 
Each  slice  was  polished  with  sand  paper.  The  CP5-4290  material  was 
calibrated  for  its  fringe  constant  by  means  of  a beam  loaded  in  pure 
bending  using  dead  weights  and  subject  to  the  same  thermal  cycle  as  for 
the  cylinders. 
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NOTE:  TRANSVERSELY— FLAWEO  CYLINDER  SUBJECTED 
SEPARATELY  TO  BOTH  UNIAXIAL  EXTENSION  AND 
INTERNAL  PRESSURE  LOADING.  THE  CYLINDERS 
WITH  INTERNAL  PRESSURE  LOADING  WERE 
SUBJECT  TO  FIXED  END-CAP  CONDITIONS. 

Figure  2.  Transverse  (A)  and  longitudinal  (B)  flaw  loading 
geometry  for  a hollow  cylinder. 

To  Improve  resolution  for  analysis,  slices  were  placed  In  an  oil 
bath  consisting  of  71.1%  by  volume  of  Halowax  oil  and  24.S%  mineral 
oil.  The  index  of  refraction  of  the  oil  was  the  same  as  the  CP5-4290 
material;  therefore,  light  scatter  was  minimized.  The  slices  were 
observed  in  a comparator  polariscope  at  10X  magnification.  By  means 
of  an  x-y  table  on  the  polariscope,  points  on  the  slices  could  be 
located  to  within  ±0.0001  in.  Fractional  fringe  orders  were  obtained 
using  Tardy  compensation.  Figure  5 shows  a photograph  of  a typical 
fringe  pattern  at  a slice  crack  tip. 


— -I 

I 

(a)  TRANSVERSE  FLAW 


Figure  S.  Typical  fringe  pattern  of  a slice. 


IV.  RESULTS  AND  DISCUSSION 

Specimen  dimensions  and  test  parameters  for  each  of  the  three 
series  of  test  are  shown  in  Tables  l,  2,  and  1.  The  notation  for  the 
test  parameters  given  in  the  tables  is  defined  in  Figure  h.  The  slices 
were  examined  in  the  comparator  polariscope.  Fringe  orders  were  taken 
at  intervals  along  a line  through  the  crack  tip  normal  to  the  crack  at 
1J1  = n / 2 . 

F.quation  (b)  must  he  used  to  resolve  the  data  for  determining  . 

This  equation  will  be  valid  only  in  a given  region  near  the  crack  tip. 
The  stress  at  points  too  near  the  crack  will  be  disturbed  by  the  finite 
crack  width  and  large  crack  tip  displacements  while  at  points  remote 
from  the  crack  tip,  the  singular  st resses  represented  bv  K(  will  fail 

to  dominate  the  total  stress  field.  Thus  a region  of  validity  in  which 
F.quation  (6)  is  assumed  to  apply  must  he  determined.  This  was  done  |1| 
by  placing  an  upper  and  lower  limit  on  the  distance  r bevond  which  no 
data  points  were  used. 
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TABLE  2.  SPECIMEN  DIMENSIONS  AND  TEST  PARAMETERS  FOR  LONGITUDINAL- FLAWED  CYLINDER 
LOADED  WITH  INTERNAL  PRESSURE 


Vfringe  / 
Crack  Width  (in 


TABLE  3.  SPECIMEN  DIMENSIONS  AND  TEST  PARAMETERS  FOR  TRANSVERSELY-FLAWED  CYLINDER 
LOADED  WITH  INTERNAL  PRESSURE 


INSIDE  SURFACE 


Figure  6.  Notation  for  the  part-circular  surface  flaw. 


The  approach  here  will  depart  from  that  in  [1],  in  favor  of  the 


method  used  by  Smith  et  al.  [5].  Solving  Equation  (6)  for  \ and 

max 

truncating  the  results  to  the  same  order  in  r as  Equation  (2)  leads  to 
the  simple  result: 


T 

max 


(7) 
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Define  an  apparent  value  of  K,  K , as,  K ■ /8wr  i . Then 

3p  1H3X 

Equation  (7)  can  be  written  as 

K * KT  + /8irr'  o 
ap  l on 


(8) 


This  equation  shows  that  in  the  region  of  valid  data,  i.e.,  the  region 
dominated  by  the  singular  stresses,  there  is  a linear  relationship 
between  the  apparent  K and  the  square  root  of  r.  Departure  from  linearity 
indicates  data  outside  the  region  of  validity. 


In  determining 
1/2 


for  a given  slice  the  values  of  K are  plotted 

versus  r . Data  points  approximately  yielding  a straight  line  are 
selected  and  all  others  are  rejected.  A least-squares  curve  is  then 
fit  to  the  selected  points  and  the  value  of  K^.  is  determined  by  taking 


K at  r = 0,  since  K = K,  at  r = 0,  Figure  7 gives  an  example  for  a 
ap  ap  I 

typical  slice.  When  the  same  points  are  used,  this  method  yields  the 
same  results  as  the  statistically  more  complex  method  of  [1].  The 
advantage  of  the  present  method  is  its  simplicity  and  the  ease  of 
selecting  data  representative  of  the  singular  stresses.  Appendix  R 
contains  a sample  listing  of  the  program  used  to  analyze  photoelasticity 
data  on  a PDP  11/40  minicomputer  with  an  RT11  software  system.  Tables  1, 
2,  and  3 and  Figures  8 through  13  summarize  the  results  of  the  testing 
program. 


Figures  8 and  9 give  the  results  for  circumferential  flaws  in  exten- 
sion for  R = 0.875  in.  and  1.500  in.,  respectively.  As  the  A/T  aspect 
ratio  increases,  an  increase  in  the  nondimensional  factor  for  all  R 

is  noted.  Figures  10  and  11  give  the  same  general  results  except  that 
the  data  appear  to  have  greater  experimental  scatter.  The  results  of 
Test  2PL  are  not  shown  in  Figure  10  due  to  excessive  data  scatter. 

Figures  12  and  13  show  similar  results  for  circumferential  flawed 
cylinders  with  internal  pressure.  Test  1PC  data  are  not  presented  due 
to  cylinder  failure.  The  results  are  generally  good  for  Tests  2PC 
through  5PC;  however.  Tests  6PC  through  8PC  show  a great  deal  of  data 
scatter. 


V.  SUMMARY  AND  CONCLUSIONS 

Photoelasticity  tests  were  conducted  on  24  surface-flawed 
hollow  cylinder  specimens  loaded  in  extension  and  internal  pressurization. 
Results  from  these  tests  are  presented  for  transverse  and  longitudinal 
semi-circular  flawed  cylinders  having  two  different  flaw  radii  and 
several  A/T  penetration  ratios.  The  tests  indicate  that  for  a fixed 
flaw  radii  the  nondimensional  stress  intensity  increases  with  an 
increase  in  the  A/T  ratio.  Results  for  these  tests  coincide  in  general 
with  those  of  [1], 
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NOMENCLATURE 

© TEST  DATA  USED  TO 
DETERMINE  K, 


© TEST  DATA  REJECTED 


Figure  11.  Stress  intensity  factor  versus  0 for  the  longitudinal 
flawed  cylinder  loaded  with  internal  pressure  (R  » 1.50  in.). 
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CIRCUMFERENTIAL  FLAW  WITH  INTERNAL  PRESSURE 
R - 0.875  In. 
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Figure  12.  Stress  intensity  factor  versus  0 for  the  transversely- 
flawed  cylinder  loaded  with  internal  pressure  (R  = 0.875  in.). 
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Appendix  A.  TEST  DATA 


This  appendix  contains  a summary  of  all  the  test  data  presented 
in  this  report.  The  test  nomenclature  is  as  follows: 

T = Extension  loading  test,  transverse  flaw. 

PL  = Internal  pressure  loading  test,  longitudinal  flaw 

PC  = Internal  pressure  loading  test,  circumferential  flaw 
(transverse  flaw) 

The  test  nomenclature  follows  the  test  number.  For  example,  1PL  refers 
to  the  first  test  run  of  a flaw  in  the  longitudinal  orientation  sub- 
jected to  internal  pressurization. 

Test  No.  IT 


Test  No.  4T 


A/T  - 0.674 

R = O.i 

875  in. 

o *9.84  psi 
m 

Slice  No. 

Angle 

(deg) 

KI 

Ki 

2/tt  o (/Et) 
m 

1 

0 

12.0330 

1.158 

3 

20 

14.8617 

1.431 

4 

20 

13.7995 

1.329 

5 

30 

13.1259 

1.264 

Test  No.  5T 


A/T  - 0.201 

R - 1. 

50  in. 

a = 9.48  psi 
m 

Slice  No. 

Angle 

(deg) 

KI 

Ki 

2/tt  a (/ttR) 
m 

1 

0 

5.4174 

0.413 

2 

10 

5.1066 

0.389 

3 

20 

4.0967 

0.312 

A/T  - 1.10 


R = 1.50  in 


a ■ 5.162  psi 
m r 


Slice  No. 

Angle 

(deg) 

KI 

KI 

2/tt  a (/ItR) 
m 

1 

15 

16.5615 

2.321 

2 

20 

13.6390 

1.912 

3 

30 

12.061A 

1.691 

A 

AO 

11.7697 

1.6A9 

Test  No.  1PL 


A/T  « 0.209 

R = 0. 

875  in. 

o =*  13.138  psi 
m 

Slice  No. 

m 

KI 

Ki 

2/tt  a (/ttR) 
m 

1 

0 

9.7267 

0.701 

2 

10 

8.890A 

0.6A1 

3 

15 

10.3859 

0.7A9 

A 

20 

6.0037 

O.A33 

5 

25 

6.9A83 

0.501 

31 


J 


Test  No.  2PL 


A/T  - 0.388 

R = 0. 

875  in. 

o * 11.898  psi 
m 

Slice  No. 

a 

K1 

Ki 

2/tt  o ( ^ttR) 
m 

1 

0 

11.0324 

0.878 

2 

10 

9.7557 

0.777 

3 

15 

14.7180 

1.172 

4 

20 

12.8083 

1.020 

5 

30 

7.7985 

0.621 

6 

35 

7.6696 

0.611 

Test  No.  3PI. 


A/T  = 0.553  R = 0.875  in.  o = 10.585  psi 

m 


Slice  No. 

Angle 

(deg) 

K1 

KI 

2 / it  o (*7r) 
m 

1 

0 

14.1988 

1.271 

2 

15 

14.9406 

1.337  j 

3 

20 

14.1567 

1.267 

4 

30 

13.8462 

1.239 

5 

35 

12.2796 

1.099 

Test  No.  4PL 


A/T  - 0.674 

R = 0, 

.875  in. 

o = 9.628  psi 
m r 

Slice  No. 

Angle 

(deg) 

KI 

Ki 

I 2/tt  a (»/ttR) 

m 

1 

0 

14.4512 

1.422 

2 

15 

13.3491 

1.314 

3 

20 

14.3452 

1.412 

4 

30 

15.2376 

1.499 

5 

35 

15.0910 

1.485 

Test  No.  5PL 


A/T  = 0.202 

R = 1. 

50  in. 

a = 10.474  psi 
m 

Slice  No. 

Angle 

(deg) 

KI 

Ki 

2/tt  a (^tR) 
m 

1 

0 

8.4531 

0.584 

2 

5 

10.0758 

0.696 

3 

10 

8.3307 

0.576 

4 

15 

7.7810 

0.538 

5 

20 

5.5351 

0.382 

33 
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Test  No.  6PL 


A/T  = 0.448 

R - 1. 

5 in. 

a = 9.494  psi 
m r 

Slice  No. 

Angle 

(deg) 

*1 

Ki 

2/ir  a ( /ttR) 
m 

1 

0 

10.8909 

0.830 

2 

10 

10.9049 

0.831 

3 

15 

10.4457 

0.796 

4 

20 

10.1675 

0.775 

5 

25 

8.5832 

0.654 

6 

30 

7.2234 

0.551 

Test  No.  7 PL 


A/T  = 0.580  R = 1.5  in.  o = 8.06  psi 

m 


Slice  No. 

Angle 

(deg) 

K 

I 

KI 

2/tt  a (/^R) 
m 

1 

0 

14.0269 

1.259 

2 

10 

14.2013 

1.275 

3 

15 

12.7411 

1.144 

4 

20 

15.5208 

1.393 

5 

25 

17 . 2004 

1.544 

7 

35 

14.2713 

1.281 

34 


— _ — - m-mr TLI 


0 

4.7515 

0.8546 

15 

5.0981 

0.9169 

20 

5.4732 

0.9844 

25 

5.0620 

0.9104 

30 

5.6509 

1.0163 

38 

5.8497 

1.0521 

Test  No.  4PC 


A/T  = 0.674  R = 0.875  in.  a = 9.91  psi 


Angle 

(deg) 

KI 

KI 

2/tt  o (/ttR) 
m 

0 

10.3759 

0.9920 

10 

11.0626 

1.0576 

15 

11.1789 

1.0687 

20 

12.5803 

1.2027 

Test  No.  5PC 


A/T  = 0.200 

R = 1 

500  in. 

o = 11.13  psi 
m 

Slice  No. 

m 

KI 

Ki 

2/tt  a (/ttR) 
m 

1 

0 

9.0737 

0.5896 

2 

10 

8.6878 

0.5645 

3 

15 

9.1832 

0.5967 

4 

20 

9.0818 

0.5901 

Test  No.  6PC 


A/T  = 0.443 

R = 1 

.500  in. 

o = 9.92  psi 
m 

Slice  No. 

Angle 

(deg) 

KI 

Ki 

2/tt  o ( /ttR) 
m 

1 

0 

8.4709 

0.6179 

2 

10 

8.7018 

0.6347 

3 

15 

9.1253 

0.6656 

4 

20 

7.0955 

0.5175 

37 


Test  No.  7 PC 


A/T  - 0.684 

R ■ 1. 

.500  in. 

a * 8.68  psi 
m 

Slice  No. 

Angle 

(deg) 

KI 

Ki 

2/ir  a (/irR) 

ID 

1 

0 

8.9268 

0.7445 

2 

10 

10.6402 

0.8874 

3 

15 

9.6247 

0.8027 

4 

25 

9.2004 

0.7673 

5 

30 

11.8390 

0.9874 

Test  No.  8PC 


A/T  = 0.940 

R = 1 

500  in. 

a =7.54  psi  j 

m r 

Slice  No. 

Angle 

(deg) 

KI 

Ki 

2 / t r a (/ttR) 
m 

1 

15 

10.2689 

0.9845 

i 2 

20 

10.1580 

0.9739 

3 

30 

11.8707 

1.1381 

4 

40 

10.4856 

1.0053 

38 


Appendix  B.  COMPUTER  CODE 


The  computer  code  shown  on  the  following  pages  was  used  to  analyze 
the  photoelasticity  data.  The  code  uses  the  technique  presented  in 
Section  IV  with  a least  squares  straight-line  curve  fit  of  the  experi- 
mental data. 


C PHOTOELASTICITY  COPE -CYLINDERS 

DIMENSION  AN ( 50 ) > AR  ( 50 ) ■:  AT  ( 50  > » AK  ( 50 ) » ID ( 50 ) 
WRITE (5 >22) 

22  FORMAT < ' NO.  OF  SLICES?') 

READ (5.1.  ) N 

1 FORMAT (13) 

DO  19  1=1 » N.l 
WRITE ( 5 » 2 > 

2 FORMAT ( ' NO.  OF  DATA  POINTS?') 

READ (5.1)  M 

WRITE (5.23) 

23  FORMAT ('  FRINGE  CONSTANT?') 

READ (5. 24)  F 

24  FORMAT (FI 0.0) 

WRITE (5.25) 

25  FORMAT ( ' SLICE  THICKNESS?') 

READ (5, 26)  T 

26  FORMAT ( F10 . 0 ) 

WRITE (5. 27) 

27  FORMAT ('  INPUT  R AND  N-F4 . 0 » F6 . 0 X ' ) 

DO  4 J=1 » M » 1 

READ (5.3)  AR ( J ) > AN ( J ) 

3 F0RMAT(F4,0.F6.0) 

4 CONTINUE 
AM  AX-- 0 . 

DO  5 J=1 > M » 1 

AT( J)=F*AN( J)/(2.*T) 

AK  ( J ) =~AT  ( J ) YSQRT  ( 8 . *3 . 1 4 1 59* AR  ( J ) ) 

IF ( AK ( J ) .GT.AMAX)  AMAX-AK ( J ). 

5 AR ( J ) -SORT ( AR ( J ) ) 

6 CALL  IP0KE(  *17041,0.  "1) 

CALL  IP0KE( "170410. *0) 

DO  7 J=1 » M » 1 

I Y= INT  < AK ( J ) * 1 000 . /AMAX  > 

IX=INT(AR( J)*1000,/AR(M) > 

CALL  I POKE ( * 1 704 1 2 . IX ) 

CALL  IP0KE( *170414. IY) 

7 CALL  IP0KE( *170414. *0) 

I TEST-' I PEEK  ( *177570) 

IF  ( I TEST . EG) . 0 ) GOTO  6 
WRITE (5. 8) 
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FORMAT ('  NO.  QF  DELETE  SAMPLES?') 

READ ( 5 . 9 ) ND 

9 FORMAT (13) 

IF(ND.EQ.O)  GOTO  20 
DO  11  J-l.NDfl 
READ (5.10)  I D ( J > 

10  FORMAT ( 13 ) 

11  CONTINUE 

20  CONTINUE 
X1=0. 

X2=0 . 

Y0=0 . 

Y1=0 . 

DO  14  J=1 » M » 1 
IF (ND.EQ.O)  GOTO  21 
DO  13  K=l.ND.l 

13  IF( J.EQ.ID(K) ) GOTO  14 

21  X1=X1+AR(J) 

X2=X2+AR< J)**2. 

YO=YO+AK  < J ) 

Y1=Y1+AK( J)*AR< J) 

14  CONTINUE 

AK1-(X2*Y0-Y1*X1  )/(FL.0AT<M-ND>*X2-Xl*Xl ) 
WRITE<5. 15)  I 

15  FORMAT < ' SLICE  N0.  = ',I3) 

DO  18  J~1 » M r 1 

AR< J)=AR< J)**2. 

WRITE (5. 16)  AR  ( J ) > AN  ( J > . AT  ( J > .« AK  < J ) 

16  FORMAT ( ' R-' ,F10.4r5X* ' N= ' F10 . 4 > 5X • 

1 ' TMAX- ' » F10 . 4 . 5X » ' KAPP= 'F10 . 4 ) 

18  CONTINUE 
WRITE <5 » 17)  AK1 

17  FORMAT  < ' K 1 = ' » F 1 4 ♦ 4 > 

19  CONTINUE 
STOP 
END 

* 
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